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ABSTRACT 

We  use  a  block  Lanczos  algorithm  for  computing  a  few  of  the  smallest 
eigenvalues  and  the  corresponding  eigenvectors  of  a  large  symmetric  matrix  rather 
than  computing  all  the  eigenvalue-eigenvector  pairs.  The  basic  Lanczos  algorithm 
generates  a  similar  matrix  which  is  block  tridiagonal  from  a  given  large  symmetric 
matrix.  The  size  of  the  generated  tridiagonal  matrix  depends  upon  the  number  of 
the  smallest  eigenvalues  to  be  computed.  The  result  is  savings  in  computations  and 
storage.  The  block  Lanczos  algorithm  is  well -suited  for  problems  involving  multiple 
eigenvalues. 

In  this  thesis,  we  develop  the  block  Lanczos  algorithm  to  estimate  the 
direction-of-arrival  (DOA)  of  a  point  source  based  on  the  observations  measured  at 
a  linear  array  of  sensors  and  compare  the  performance  with  that  of  a  single  vector 
Lanczos  algorithm.  The  results  of  the  computer  simulation  experiments  conducted 
with  this  method  are  presented  and  discussed. 
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I.    INTRODUCTION 

It  is  necessary  to  estimate  the  power  spectral  density  (PSD)  function  in 
several  signal  processing  problems.  Of  particular  interest  is  the  estimation  of 
spectral  lines  in  noise.  One  related  problem  is  the  direction— of— arrival  estimation  of 
point  sources  based  on  array  measurements.  Methods  based  on  subspace  estimation 
have  been  proposed  in  this  regard  [Ref.  9].  Much  of  the  basic  literature  on  these 
methods  has  been  borrowed  from  functional  and  numerical  analysis  [Refs.  8,  21]. 
Determination  of  the  subspaces  needs  a  complete  eigendecomposition  of  the  given 
data  autocorrelation  matrix.  The  general  eigendecomposition  of  an  nxn  matrix 
requires  0(n3)  computations.  Often  we  only  need  to  compute  a  few  of  either  the 
smallest  or  the  largest  eigenvalues  and  the  corresponding  eigenvectors  of  a  large 
symmetric  autocorrelation  matrix  rather  than  all  the  eigenpairs  of  the  matrix.  In 
this  research,  we  develop  an  algorithm  for  computing  a  few  of  the  smallest 
eigenvalues  and  the  corresponding  eigenvectors  and  apply  this  to  high  resolution 
spectral  estimations  problems. 

The  algorithm  we  develop  is  an  extension  of  the  method  of  minimized 
iterations  due  to  Lanczos  [Ref.  9].  Paige  [Ref.  4]  experimented  with  the  Lanczos 
algorithm  and  found  that  a  few  of  the  extreme  eigenvalues  of  a  tridiagonal  matrix 
would  often  converge  rapidly  to  the  similar  eigenvalues  of  a  real  symmetric  matrix 
R  much  before  the  entire  set  of  eigenvalues  are  computed. 

The  block  Lanczos  algorithm  is  an  extension  of  the  basic  single  vector  Lanczos 
algorithm.  We  iterate  with  a  block  of  vectors  rather  than  with  a  single  vector,  and 
generate   a   reduced   block    tridiagonal    matrix    that    is   similar   to    the   original 


autocorrelation  matrix  R.  The  basic  idea  is  that  the  eigenvalues  of  the  block 
tridiagonal  matrix  are  approximately  the  same  as  the  extreme  eigenvalues  of  R. 
Several  researchers  have  reported  the  block  Lanczos  algorithm  and  its  variants,  in 
particular,  Cullum  and  Willoughby  [Ref.  1],  Kahan  and  Parlett  [Ref.  6]  and  Golub 
and  Underwood  [Ref.  4]. 

The  objective  of  this  thesis  is  to  develop  a  block  Lanczos  algorithm  to  extract 
a  few  of  the  smallest  eigenvalues  and  then  estimate  the  corresponding  eigenvectors. 
The  smallest  eigenvalues  are  said  to  correspond  to  the  noise  subspace  of  the  spectral 
or  array  measurements.  The  proposed  algorithm  will  be  used  to  estimate  the 
spectral  lines  which  may  represent  the  direction-of-arrival  of  point  sources  in  low 
signal-to-noise  ratio  environments.  The  block  Lanczos  algorithm  will  have  to  be 
compared  with  the  single  vector  case  with  respect  to  the  spectral  line  estimation 
performance.  With  these  goals  in  mind,  we  now  proceed  to  discuss  how  the  thesis 
is  organized. 

In  Chapter  II  we  introduce  and  summarize  the  single  vector  Lanczos 
algorithm.  Complete  and  selective  reorthogonalization  of  Lanczos  vectors  is 
discussed  there.  In  Chapter  III,  we  describe  and  develop  the  block  Lanczos 
algorithm  and  present  the  results  of  the  experiments  carried  out  with  a  computer 
program  implementing  this  algorithm.  Also,  in  this  chapter,  we  compare  the 
performance  of  DOA  estimation  of  the  block  Lanczos  algorithm  with  that  of  the 
single  vector  Lanczos  algorithm.  Finally,  in  the  last  chapter,  we  discuss  and 
summarize  the  results  of  the  Lanczos  method  and  also  make  some  recommendations 
for  future  work. 


II.    SINGLE  VECTOR  LANCZOS  ALGORITHM 

The  single  vector  Lanczos  algorithm  is  used  to  tridiagonalize  a  real  symmetric 
matrix  R.  This  algorithm  is  based  on  the  concepts  such  as  Krylov  sequences  and 
subspaces,  orthogonal  projections  of  matrices,  and  Ritz  vectors  [Ref.  1].  In  this 
chapter  we  describe  the  general  Lanczos  recursion  and  a  practical  Lanczos 
algorithm.  Issues  related  to  the  reorthogonalization  of  Lanczos  vectors  are  also 
discussed. 

A.      LANCZOS  RECURSION 

Given  an  n*n  real  symmetrix  matrix  R  and  an  arbitrary  n*l  unit  norm  vector 
kj,  we  can  obtain  a  sequence  of  vectors  defining  the  n  dimensional  Krylov  subspace 
as  follows  [Ref.  1] 

k2  =  Rk: 

k3  =  Rk2  =  R\ 


kn  =  Rk^  =  R»-%.  (2-1) 


We  can  now  form  the  Krylov  matrix  of  rank  m  as  follows 


^m  —  1*1   k2  k3    •    •    •    kj 

=  [kl5  Rkl5  R2k1:  •  •  •,R»-1kJ  (2-2) 


and  the  Krylov  subspaces  Xm,  m=l,2,-  •  -,n,  are  then  defined  as 

JTm(R,k,)  =  span  {lq,  Rkl3  R2kl5  •  •  -,  R^kJ.  (2-3) 

We  now  proceed  to  obtain  a  tridiagonal  matrix  which  has  the  same  extreme 
eigenvalues  as  the  given  symmetric  matrix  R.  If  the  columns  of  a  matrix  Qm  are 
orthonormal,  then  the  matrix 


Tm  =  Q^RQm  (2-4) 


is  an  mxm  tridiagonal  matrix  which  is  an  orthogonal  projection  of  R  onto  the  space 
spanned  by  the  columns  of  Qm.  Lanczos  [Refs.  1,  9,  14]  proposed  a  recursion  to 
generate  the  columns  of  the  matrix  Qm  =  [qt  ^  •  •  •  qj.  The  vectors  qv  q2,  •  •  •,  qm 
are  referred  to  as  Lanczos  vectors.  Thus,  for  a  given  real  symmetric  matrix  R,  the 
Lanczos  recursion  produces  a  tridiagonal  matrix  Tm  as  a  projection  of  R  onto  the 
corresponding  subspace  span{Qm},  spanned  by  the  Lanczos  vectors  generated. 
These  subspaces  are  in  fact  Krylov  subspaces.  [Ref.  1] 

The  classical  Lanczos  recursion  for  R  starts  with  an  nxl  initial  Lanczos  vector 
qt  which  is  randomly  generated  and  normalized.  Initially  we  define  /?x=  0  and  q0=  0. 
We  then  compute  av  /?i+1  and  the  Lanczos  vectors  q±  for  i=l,2,«  •  -,m  as  follows: 


Of|  =  q^Rqi  (2-5) 


A+i<li+i  =  R(li  ~  <*ft\  ~  ^i-i  (2-6) 


A+i  =  Qi+iR(li- 


(2-7) 


The  Lanczos  matrix  Tm  is  generated  by  appropriately  arranging  ai  and  /?i+1,  where 
ai  are  the  diagonal  elements  and  /?i+1  are  the  subdiagonal  elements  of  the  tridiagonal 
matrix,  so  that 


and 


Tm(i,i)  =  orj        for  l<i<m, 


(2-8) 


TJi,i+l)  =  Tm(i  +  l,i)  =  /?i+1      for   l<i<m-l, 


(2-9) 


resulting  in 


T    = 


<*1     $2 
02    a2    h 
h    «3 


'  An 

An      Qm 


(2-10) 


The  vectors  a^  and  /^qj.j  in  Eqn(2-6)  are  the  orthogonal  projections  of  vectors 
Rqj  onto  the  two  most  recently  generated  Lanczos  vectors  qi  and  qj_r  Thus, 
updated  Lanczos  vector  qi+1  is  computed  by  orthogonalizing  the  vector  Rqj  with 
respect  to  previously  computed  Lanczos  vectors  q±  and  qw.  The  classical  Lanczos 
recursion  can  be  condensed  into  matrix  notation  by  a  group  of  Lanczos  vectors  Qm 


and  the  Lanczos  matrix  T  .   We  obtain  the  following  matrix  equation 


RQm  =   QmTm  +  /?i+1qi+,eT  (2-11) 


where  e^  is  the  unit  vector  whose  2th  element  is  1  and  whose  other  elements  are  0. 
Note  that,  in  this  recursion,  R  is  never  modified.  Also,  storage  is  needed  only  for 
the  Lanczos  vectors  qj.1?  q1  and  qi+1,  the  Lanczos  matrix  Tm,  and  the  given  matrix 
R. 

The  classical  Lanczos  procedure  attempts  to  maintain  the  orthogonality  of  the 
Lanczos  vectors.  Due  to  the  roundoff  and  other  numerical  errors,  however, 
orthogonality  between  the  Lanczos  vectors  can  only  be  maintained  by  incorporating 
some  kind  of  explicit  reorthogonalization  as  the  Lanczos  vectors  are  computed. 
Note  that  the  reorthogonalization  of  these  vectors  requires  extra  storage  for  keeping 
all  of  the  Lanczos  vectors  as  well  as  additional  computations.  We  discuss  the 
reorthogonalization  in  a  later  section. 

B.   PRACTICAL  LANCZOS  ALGORITHM 

In  this  section  we  focus  on  a  Lanczos  algorithm  with  no  explicit 
reorthogonalization  incorporated.  The  loss  of  orthogonality,  if  it  occurs,  may  result 
in  spurious  eigenvalues.  Although  the  orthogonality  of  the  Lanczos  vectors  is  lost, 
some  of  the  eigenvalues  of  R  will  siill  appear  as  eigenvalues  of  the  Lanczos  matrix  if 
we  make  the  tridiagonal  matrix  large  enough.  Since  we  are  generally  interested  in 
only  a  few  of  the  smallest  eigenvalues  and  their  corresponding  eigenvectors,  the  loss 
of  orthogonality  does  not  critically  affect  finding  approximate  eigenvalues  and 
corresponding  eigenvectors  of  R  [Ref.  1].   Nevertheless,  it  is  possible  to  find  accurate 


eigenvalues  and  eigenvectors  by  using  the  method  in  an  iterative  scheme  using  no 
reorthogonalization  at  all,  even  in  the  face  of  total  loss  of  orthogonality  [Ref.  12]. 
Now,   we  present  the  practical  single  vector  Lanczos  algorithm. 

To  generate  the  tridiagonal  Lanczos  matrix  Tm  we  should  compute  q{  and  /?i+1 
which  are  the  diagonal  and  sub-diagonal  elements  of  the  tridiagonal  matrix,  where 
i=l,2, •  •  -,m.  Given  an  n*n  real  symmetric  matrix  R,  we  start  with  an  n*l  initial 
arbitrary  vector  qj  which  is  normalized  such  that  |q1|2=l  as  in  Eqn(2— 1).  We  now 
define  an  intermediate  vector 

Ul  =  RQl  (2-12) 

and  initialize  a0=0  and  /?0=0.  The  single  vector  recursion  is  then  carried  out  for 
i=l,2,«  •  -,m  and  can  be  summarized  as  follows: 


a,  =  q]^  (2-13) 


Wj  =  Uj  -  a^  (2-14) 


Vi  =  +  /wfwi  (2-15) 


<li+i  =  wi/ri+i  (2~16) 


Pi  =  qlRqi+i  (2-17) 


ui+1  -  Rqi+1  -  ftqi  (2-18) 


where  u5  and  W;  are  the  orthonormal  projection  vectors  for  the  Lanczos  vector  qv 
Here  Eqn(2— 13)  and  Eqn(2— 17)  use  the  modified  Gram-Schmidt  orthogonalization 
to  compute  the  coefficients  a-^  and  /?i+1,  respectively.  If  qj  is  orthogonal  to  q^.v 
then  qi+1  will  theoretically  be  orthogonal  to  both  qi.l  and  qj.  This  algorithm  is 
quite  popular  despite  the  requirement  for  small  amounts  of  extra  computations 
[Ref.  12].  Now  the  tridiagonal  matrix  Tm  is  generated  by  simply  filling  it  with  ai 
and  /?j  for  its  entries  as  shown  in  Eqn(2— 10). 

1.       Eigenvalue  Computation 

To  find  the  eigenvalues  /ij  of  the  mxm  Lanczos  matrix  Tm,  we  may  use 
the  bisection  method  and  Sturm  sequencing  [Ref.  5].  Actually,  one  could  obtain 
both  eigenvalues  and  eigenvectors  of  Tm  by  using  such  methods  as  the  QR  algorithm 
[Ref.  16].  However,  since  we  need  only  a  few  of  the  smallest  eigenvalues  of  Tm,  we 
choose  the  bisection  method  to  compute  them  as  detailed  in  the  following. 

Given  the  tridiagonal  matrix  Tm,  we  define  the  characteristic 
polynomials  p0(//),  p^/x),  •  •  •,  pm(//)  as 

PoM  =  1  (2-19) 

Pj(//)  =  det(Tm-//I),  (2-20) 

for  j=l,2,«*',m.  Expanding  the  determinant  in  Eqn(2— 20)  yields  the  recursive 
expression  [Ref.  5:  pp.  305-307]: 

Pj(/i)  =  (ftj  -tiVyM  -/3fPj-2M,       j=2,3r  •  -,m.  (2-21) 


The  zeros  of  the  polynomial  pm(//)  are  the  eigenvalues  of  the  tridiagonal  matrix  Tm. 
Here  we  are  interested  in  only  a  few  of  the  smallest  eigenvalues  of  R.  In  order  to 
compute  these  values  we  first  define  a  range  [x,  y]  in  which  all  the  desired 
eigenvalues  lie.  We  then  carry'  out  the  following  iteration  to  implement  the  Sturm 
sequencing  property.   If  p(x)p(y)  <  0  and  x  <  y,  then  the  iteration 

|x-y|  >  c(  |x|  +  |y|  ) 


u-  x  +  y 


y  =  i*    if  pmWpm(y)  <  o  (2-22) 

x  =  H     if  pm(x)pm(y)  >  0 

is  guaranteed  to  converge  to  a  zero  of  pm(//),  i.e.,  to  an  eigenvalue  of  Tm.  The  value 
i  is  the  machine  unit  roundoff  error  and  the  limits  on  the  range  [x,  y]  are  given  by 

x  =  min  ftj  —  |  /?j  |  —  |  /?}_!  | 

(2-23) 
y=  max  a{  +  \/3{\  +  l/VJ 

i 

where  we  have  /30  —  /?m+1  =  0.  The  bisection  method  computes  the  eigenvalues  with 
small  relative  error,  regardless  of  their  magnitude.  This  is  in  contrast  to  the 
tridiagonal  QR  iteration,  where  the  computed  eigenvalues  can  have  only  small 
absolute  error. 


2.       Eigenvector  Computation 

Now,  we  need  to  compute  the  eigenvectors  of  R  to  complete  the 
eigenpair  estimation  problem.  There  are  two  techniques  to  compute  the 
approximate  eigenvectors  of  R.  The  first  technique  uses  the  Ritz  vector  X:  and  the 
corresponding  algorithm  uses  the  relationship  [Ref.  14] 

xj  =  Qmzj,       j=l,2,...,s  (2-24) 

where  z-  is  one  of  the  eigenvectors  of  Tm,   Qm  is  the  matrix  of  Lanczos  vectors,   and 

The  second  technique  to  obtain  the  eigenvectors  of  R  corresponding  to 
the  selected  eigenvalues  of  Tm  uses  the  Rayleigh  quotient  iteration.  The  iteration 
makes  the  Rayleigh  quotient  r(xk)  converge  to  /Zj  which  is  one  of  the  smallest 
eigenvalues  of  Tm.  The  Rayleigh  quotient  of  an  eigenvector  xk  which  minimizes 
|(R-rkI)xk||2isgivenby 


xkRxk 


rk  =  r(xk)  = = .  (2-25) 


xkxk 


As  r(xk)  approaches  an  eigenvalue  /^  of  Tm,  then  the  solution  to  (R  —  rkI)xk  =  bk 
will  be  an  approximate  eigenvector  by  using  the  inverse  iteration  theory  where  bk  is 
a  vector  close  to  zero  [Ref.  14]. 

We  now  briefly  present  the  inverse  iteration  algorithm.    First,    we  pick 
an  arbitrary  unit  vector  Xq;   then  the  iteration  proceeds  as  follows 


10 


fork  =  0,1, 


r^=x^TIbc^,  x^Tx^=l 


solve  (R  -  rkiJ  I)yk+1  =  xk!)    for  yk+1 


(2-26) 


xk+i  —  yk+i/Byk+iB2- 

If  rkx)  converges  to  one  of  the  smallest  eigenvalues  n{  of  Tm  where  i=l,2,«  •  -,s,  then 
we  stop  the  iteration  and  the  corresponding  vector  x{kl\  is  the  required  eigenvector 
of  R. 

Having  discussed  the  Lanczos  algorithm  and  the  methods  to  compute  the 
desired  eigenpairs,  we  now  proceed  to  address  some  issues  related  to  the 
reorthogonalization  of  the  Lanczos  vectors. 

C.      REORTHOGONALIZATION 

As  mentioned  in  the  previous  sections,   the  Lanczos  vectors  ql5  q2,  •  •  • ,  qm  lose 

mutual  orthogonality  as  the  number  of  steps  m  in  the  algorithm  increases.     The 

T 

requirement  that  QmQm  =  Im  is  then  destroyed  by  the  roundoff  errors  and  the 

algorithm  is  described  as  unstable.  A  few  steps  later,  the  matrix  of  Lanczos  vectors 
Qm  may  not  even  have  full  rank,  i.e.,  the  Lanczos  vectors  may  become  linearly 
dependent.  As  a  result  there  is  no  guarantee  that  Tm  will  bear  any  useful 
relationship  to  R.  And  the  orthogonality  among  vectors  ql5  q2,  •  •  • ,  qm  disappears. 
The  ideal  Lanczos  algorithm  should  terminate  (fim  =  0)  for  some  m  <  n,  but  in 
practice  the  process  goes  on  forever  computing  more  and  more  spurious  eigenvectors 
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for  each  correct  eigenpair  it  discovers.  [Ref.  14] 

Thus,  sometimes  we  need  to  maintain  the  orthogonality  of  the  Lanczos  vectors 
for  finding  more  accurate  eigenvalues  and  eigenvectors  of  R  and  to  avoid  the 
computation  of  any  spurious  pairs.  We  now  add  the  reorthogonalization  step  in  the 
single  vector  Lanczos  procedure  presented  in  Section  A  and  B  of  this  Chater.  There 
are  two  techiques  for  reorthogonalization,  namely,  complete  reorthogonalization 
and  selective  reorthogonalization.  These  are  discussed  in  the  following. 

1.       Complete  Reorthogonalization 

Complete  reorthogonalization  incorporates  a  Householder  matrix 
computation  into  the  Lanczos  algorithm  for  producing  Lanczos  vectors  that  are 
orthogonal  to  within  the  working  accuracy.  This  is  effective  at  maintaining  the 
stability  of  the  system.  The  following  step  is  inserted  into  the  Lanczos  algorithm 
(Eqns(2— 13)— (2—18))  after  computing  a  projection  vector  Wj: 


wj  =  wj  -  qj(qjWi),  forj=i,i-l,--.,2,l  (2-27) 


thus  Wj  is  explicity  orthogonalized  against  q}  and  qj.j.  If  a  Householder  matrix  Pj  is 
determined  so  that  [P0  Pj  •••  Pj]  [w0,  wl5  •••,  wj  is  upper  triangular,  then  it 
follows  that  the  (t+l)st  column  of  [P0  •••  Pj  is  the  desired  unit  vector.  An 
example  of  a  complete  reorthogonalization  Lanczos  scheme  is  summarized  below. 

T      T 

First,    we  determine  the  initial  Householder  matrix  P0  =  I  —  v0v0/v0v0    so  that 

T 

P0w0  =  er    From  Eqn(2— 5)  we  can  then  obtain  e^  =  q  j  Rqj  and  implement  the 

following  recursion  for  i=l,2, *  •  -,m— 1:  [Ref.  5:  pp.  334—345] 
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Fi  =  (R  -  ajltai  -  ft^qj.,     (/?0q0=0) 


w  =  (Pi-!  •  • '  PoK  (2-28) 


T       T 

determine  P1  =  I  —  2VJVJ/VJVJ 


p.Wj  =  (Wl,  •••,  wi5  /?j,  0,  •••,  0)T 


Qi+i  =  (po  •  •  *  pi)ei+i 


»i*i  =  qT+iRQi+i- 


In  Eqn(2— 28)  the  computed  Lanczos  vectors  q{  are  now  mutually  orthogonal  to  the 
working  precision  which  follows  from  the  roundoff  properties  of  Householder 
matrices  [Ref.  5].  However,  the  complete  reorthogonalization  of  Lanczos  vectors 
requires  extra  computations  and  all  the  computed  vectors  q{  q2  •  •  •  q^  wt  w2  •  •  ■  Wj 
need  to  be  stored.   It  negates  any  advantage  of  the  Lanczos  algorithm. 

2.       Selective  Reorthogonalization 

Selective  reorthogonalization  is  computationally  more  efficient.  A  small 
modification  to  the  exact  version  of  the  Lanczos  algorithm  ensures  that  the  Lanczos 
vectors  qv  q2,  q3'--  maintain  the  orthogonality.  At  the  ith  step  the  Lanczos 
algorithm  (Eqns(2— 13)  —  (2—18))  produces  the  matrix  of  Lanczos  vectors  Qv  the 
tridiagonal  matrix  Tj=  QjRQ^  and  the  residual  vector  Wj  =  (RQj  —  QjT^ej.  The 
equation  to  compute  the  residual  vectors  RQj  =  Q{^\  +  w^  can  then  be  rewritten 
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as 


RQi  =  QiTj  +  wieT  -  Fi5  (2-29) 


where  Fi  accounts  for  the  roundoff  errors.    When  the  Lanczos  vectors  lose  their 
orthogonality  we  have 


|I-QlQil<*i  (2-30) 


where  it  is  required  to  keep  k±  <  k  for  some  k  in  the  interval  (n6,  0.01)  and  t  is  the 
relative  precision  of  the  arithmetic. 

The  loss  of  orthogonality  goes  hand  in  hand  with  the  convergence  of  a 
Ritz  pair.  Suppose  that  the  symmetric  QR  algorithm  [Ref.  16]  is  applied  to  Tj  and 
renders  the  computed  Ritz  values  dv  •  •  • ,  61  and  a  nearly  orthogonal  matrix  of 
eigenvectors  Zj=[zj,z2,-  •  -,zj  of  Tj.  Then  the  Ritz  vectors  which  are  the 
approximate  eigenvectors  of  R  are  given  by 

Xj  =  1x^X2,  ...,xi]  =  QiZi.  (2-31) 

We  can  show  that  the  absolute  value  of  the  inner  product  between  qi+1  and  X:  is 


T  t  |R|2 

1^' g   |/),|  L.,1  (2"32) 


where  the  denominator  term  can  be  approximated  as 


l^llzi^lRxj-^XjIa  (2-33) 
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for  j=l,  2,  •--,  i.  The  recently  updated  Lanczos  vector  qi+1  is  forced  to  have  a 
unwanted  component  in  the  direction  of  any  converged  Ritz  vector.  Consequently, 
instead  of  orthogonalizing  qi+1  against  all  of  the  previous  Lanczos  vectors,  we  can 
achieve  the  same  effect  by  orthogonalizing  it  against  the  converged  Ritz  vectors 
[Ref.  5].  A  selective  reorthogonalization  method  based  on  this  technique  is 
discussed  in  Parlett  and  Scott  [Ref.  15].  In  this  method,  a  computed  Ritz  pair  (#,x) 
is  considered  a  good  approximation  if  it  satisfies 

lR*-0xl2:^lRfl2  (2-34) 

where  e  is  the  machine  precision  constant.  After  computing  qi+1,  it  is 
orthogonalized  against  each  good  Ritz  vector  [Ref.  5].  Selective  reorthogonalization 
prevents  the  computation  of  many  spurious  eigenvectors.  This  means  that  the  extra 
computations  and  storage  required  in  this  method  are  less  than  those  required  to 
implement  the  complete  reorthogonalization  since  there  are  usually  fewer  Ritz 
vectors  than  Lanczos  vectors.  When  all  of  the  eigenvalues  of  R  are  required, 
however,  then  the  selective  reorthogonalization  procedures  are  too  expensive  to 
implement  [Ref.  12]. 

D.      RESULTS 

In  this  section  we  will  discuss  the  experimental  results  based  on  computer 
simulation  of  the  single  vector  Lanczos  algorithm  for  estimating  the 
direction-of-arrival  of  point  sources. 
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1.       Experimental  Set  Up 

Consider  that  we  receive  signals  at  a  linear  array  containing  M  equally- 
spaced  sensors.  The  signal  is  modeled  as  a  sum  of  /  sinusoids,  the  normalizing 
spatial  frequencies  of  which  are  propotional  to  their  bearings  from  0.0  (=  0°)  to 
0.5  (=  90°),  in  additive  random  noise  with  fixed  variance.  The  received  signal  at 
the  mth  sensor  location  is  then  given  by 

/ 
z(m)  =  ^  A1cos{27rfim)  +  n(m)  (2-35) 

i  =  l 

where  Ai  is  the  amplitude  of  the  ith  sinusoid,  /j  is  the  normalized  spatial  frequency 
of  the  zth  signal  which  represents  the  bearing,  and  n(m)  is  zero-mean  white 
Gaussian  noise  with  variance  a2.  The  amplitude  of  the  signal  A^  and  the  noise 
variance  o2  will  determine  the  signal  to  noise  ratio  of  the  Ith  signal, 


SNRj  =  10  log( ).  (2-36) 

a2 


Figure  1  shows  the  block  diagram  of  the  direction -of- arrival  estimation 
algorithm  considered  in  this  work.  After  measuring  the  signals  received  at  each 
sensor,  the  sensor  output  is  passed  through  a  bank  of  bandpass  filters.  This  helps 
prefilter  the  noise  over  the  selected  frequency  band  of  frequencies  and  provides  some 
initial  processing  gain  [Ref.  22].  The  autocorrelation  matrix  is  then  computed  by 
taking  data  from  the  outputs  of  corresponding  bandpass  filters  at  each  sensor.  An 
n*n  autocorrelation  matrix  is  generated  as  follows 
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M-l-k 


Rzz(k)  =  ^r    I    z(m+k)zH  forO<k<n-l.  (2-37) 


m  =  0 


We  now  have  computed  an  n*n  autocorrelation  matrix  of  z(m),  Rzz,  by 
using  M  data  samples.  The  eigenvectors  xi  of  Rzz  corresponding  to  the  smallest 
eigenvalues  are  computed  by  using  the  single  vector  Lanczos  algorithm.  The  power 
spectral  density  estimates  are  computed  as: 


Sii}(/)  = 


1 


n-  1 


j  =  l 


z=eJ 


j27T/ 


(2-38) 


th 


where  x^  are  the  elements  of  the  fn  eigenvector  Xj 


2.       DOA  Estimation 

A  few  of  the  extreme  eigenvalues  and  their  corresponding  eigenvectors  of 
a  large  real  symmetric  matrix  R  can  be  obtained  by  using  the  single  vector  Lanczos 
algorithm.  Those  eigenpairs  are  the  approximate  eigenpairs  of  R.  Each  eigenvector 
has  spectral  information  to  determine  the  bearing  of  a  source  relative  to  an  array  of 
sensors.  The  physical  implementation  of  a  DOA  estimation  scheme  is  shown  in 
Figure  1.  The  correlator  generates  an  autocorrelation  matrix  of  the  received  signal. 
The  Lanczos  algorithm  and  the  eigendecomposition  produce  the  noise  subspace 
eigenvectors  that  estimate  the  spatial  power  spectral  density  (PSD)  function,  given 
in  Eqn(2— 38).  The  PSD  function  represents  the  direction -of-arrival  of  point 
sources  as  spectral  peaks.  In  this  thesis  we  have  used  three  different  ways  of 
estimating  the  PSD  function:     the  individual  eigenvector  spectra,     the  algebraic 
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averaging  of  a  few  eigenvectors,  and  the  spectral  multiplication  of  the  individual 
spectra. 

The  computer  simulation  experiment  consists  of  an  equally -spaced  array 
of  100  sensors  arranged  in  a  linear  fashion.  By  using  a  filter  bank  as  shown  in 
Figure  1,  the  signals  have  known  temporal  frequency  with  unknown  bearings.  The 
autocorrelation  matrix  size  is  chosen  to  be  25x25.  The  number  of  eigenpairs 
computed  is  chosen  according  to  the  number  of  iterations  in  the  Lanczos  algorithm. 
We  have  used  5  dB,  0  dB,  —5  dB  and  —10  dB  SNR  cases  for  the  direction-of-arrival 
estimation.  The  results  indicate  the  ability  of  the  algorithm  to  determine  the 
number  of  targets  and  bearing  resolution  for  various  directions  and  different  SNRs. 
In  each  case  we  have  used  the  five  smallest  eigenvalues  and  their  corresponding 
eigenvectors  for  computing  the  PSD  function.  (We  can,  however,  choose  more 
eigenvectors  at  the  expense  of  more  computations.) 

We  can  use  eigenvector  averaging,  which  is  the  algebraic  averaging  of 
the  computed  eigenvectors  corresponding  to  the  smallest  eigenvalues  of  R.  It  has 
improved  the  performance  compared  to  that  of  the  individual  eigenvector  spectra. 
Further  improvement  in  results,  however,  can  be  obtained  by  spectral  multiplication 
of  the  individual  eigenvector  spectra.   The  estimate  is  given  by 


Sxx(/)  =  rrs*j,(/)  (2-39) 

i  =  l 


where  J  is  a  predetermined  number  ( J  <  m  <  n). 

We  now  consider  several  examples  to  study  the  estimation  performance 
of  the  Lanczos  algorithm  and  the  consequent  eigenpair  computation.  Most  of  the 
results  use  the  Ritz  vector  method  with  no  reorthogonalization.  Example  1  is  the 
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detection  of  a  single  target  at  9°  for  different  SNRs.  Figure  2  shows  the  overlayed 
individual  spectra  of  the  five  eigenvectors  corresponding  to  the  smallest  eigenvalues 
at  an  SNR  of  5  dB.  Note  that  the  spectrum  of  each  eigenvector  has  several  spurious 
peaks,  but  each  eigenvector  has  a  common  peak  at  the  true  bearing  9°.  Figure  3(a) 
shows  a  plot  of  the  PSD  function  which  is  the  average  of  5  eigenvectors.  The 
averaging  improves  the  estimation  performance  considerably.  However,  further 
improvement  was  obtained  by  using  the  spectral  multiplication  as  defined  in 
Eqn(2— 39)  and  the  result  is  shown  in  Figure  3(b).  As  can  be  seen,  even  though  it 
has  two  small  spurious  peaks,  it  has  greatly  improved  the  nulls  and  the  peak  at  9°. 
In  the  remainder  of  the  thesis,  we  use  the  spectral  multiplication  method  to 
compute  the  PSD  function.  Figure  4  and  5  show  the  results  at  0  dB  and  —5  dB  SNR 
respectively.  More  spurious  peaks  are  observed  as  the  SNR  is  lowered.  At  an  SNR 
of  —10  dB  (see  Figure  6(a)),  several  spurious  peaks  are  seen  which  are  almost  as 
large  as  the  true  bearing.  Improved  performance  is  obtained  as  shown  in  Figure 
6(b)  by  using  more  eigenvectors  (7  eigenvectors)  in  this  case  in  the  spectral 
multiplication. 

In  Example  2  we  have  3  targets  at  34°,  36°  and  54°.  Notice  that  two 
targets  are  very  closly  spaced  in  bearing.  Figure  7(a)  is  obtained  by  using  the  Ritz 
vector  method  and  Figure  7(b)  is  obtained  by  using  the  Rayleigh  quotient  iteration 
to  compute  the  eigenvectors  at  0  dB.  Although  the  two  targets  are  very  closely 
spaced,  good  resolution  is  clearly  achieved  and  almost  no  spurious  peaks  are  seen  in 
both  results.  Figure  8(a)  and  Figure  8(b)  are  obtained  by  using  no 
reorthogonalization  and  with  complete  reorthogonalization  respectively  at  an  SNR 
of  —5  dB.  As  mentioned  earlier,  loss  of  orthogonality  does  not  affect  the  results  for 
a   few   eigenvectors    in    the   single   vector    Lanczos    procedure.       The    result    of 
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reorthogonalization  is  almost  the  same  as  the  result  of  no  reorthogonalization. 
Resolution  is  achieved  in  both  cases  at  —5  dB,  but  several  spurious  peaks  are 
present.  Figure  9(a)  shows  the  result  using  3  eigenvectors  while  Figure  9(b) 
indicates  the  performance  using  5  eigenvectors  at  an  SNR  of  —10  dB.  Note  that 
sufficient  spectral  resolution  is  not  acheived  to  discriminate  the  targets  located  at 
34°  and  36°.  A  number  of  spurious  peaks  are  higher  than  the  peak  at  36°  making  it 
impossible  to  accurately  determine  the  number  of  targets  as  well  as  their  locations. 

The  results  in  this  chapter  indicate  that  the  eigenvectors  found  using  the 
single  vector  Lanczos  algorithm  are  sufficiently  accurate  to  determine  the  spectrum. 
The  spectral  multiplication  scheme  achieves  the  best  spectral  estimation 
performance.  This  algorithm  provides  savings  in  computations  and  storage  because 
it  needs  to  compute  only  a  few  of  the  extreme  eigenvalues  and  eigenvectors  of  a 
large  real  symmetric  matrix.  Loss  of  orthogonality  is  not  critically  affected  when  we 
need  to  find  only  a  few  eigenvalues  and  eigenvectors  of  a  large  autocorrelation 
matrix. 
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Figure  1.       Experiment  Set  up 
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Figure  2.       One  target  at  9°  (5  dB)  :     overlay  of  5  eigenvectors 
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Figure  3.       One  target  at  9°  (5  dB)  :     (a)  average  of  5  eigenvectors 
(b)  product  of  the  PSDs  5  eigenvectors 
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Figure  4.       One  target  at  9°  (OdB)  :     product  of  the  PSDs  5  eigenvectors 


10° 

IO"1 

10  z 

to 
C 
V 

a 

io-3 

10" 

a 

*-> 

V 
IX 

m 

10"5 
10  s 

D-, 

IO"7 

10   " 

io-" 

jO-m 

40  60 

Rearing 


100 


Figure  5.       One  target  at  9°  (-5dB)  :     product  of  the  PSDs  5  eigenvectors 
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Figure  6.       One  target  at  9°  (-10  dB)  :     product  of  the  PSDs 
(a)  5  eigenvectors     (b)  7  eigenvectors 
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Figure  7.     Three  targets  at  34°,  36°  and  54°  (0  (IB)  :     product  of  the  PSDs 
5  eigenvectors  (a)  Ritz  vector  (b)  Rayleigh  quotient  iteration 
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Figure  8.     Three  targets  at  34°,  36°  and  54°  (-5  dB)  : 

(a)  no  reorthogonalizing  (b)  with  complete  reorthogonalizing 
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III.  BLOCK  LANCZOS  ALGORITHM 

In  the  previous  chapter  we  presented  the  single  vector  Lanczos  method.  The 
single  vector  method  can  be  used  to  find  a  few  extreme  eigenvalues  of  any  given  real 
symmetric  matrix  R.  However,  this  method  does  not  determine  the  multiplicities  of 
the  eigenvalues  directly.  Besides,  it  does  not  determine  a  complete  basis  for  the 
invariant  subspace  corresponding  to  any  such  multiple  eigenvalue.  Here  we  consider 
an  alternative  approach  which  directly  determines  the  multiplicities  of  the 
eigenvalues  and  the  corresponding  eigenvectors  [Ref.  1]. 

The  method  known  as  the  block  Lanczos  algorithm  is  an  extension  of  the 
Lanczos  algorithm  in  which  a  block  of  vectors  rather  than  a  single  vector  is  iterated 
[Ref.  4].  We  produce  a  block  tridiagonal  matrix  in  place  of  the  usual  tridiagonal 
matrix  produced  in  the  single  vector  Lanczos  method.  The  block  Lanczos  method 
can  be  used  in  a  manner  proposed  by  Paige  [Ref.  10].  That  is,  one  can  compute  a 
sequence  of  estimates  of  the  eigenvalues  of  the  matrix  R  from  the  block  tridiagonal 
matrix. 

Several  researchers  have  worked  on  the  block  Lanczos  algorithm,  in  particular 
Kahan  and  Parlett  [Ref.  6],  Cullum  and  Willoughby  [Ref  1],  and  Golub  and 
Underwood  [Ref.  4].  In  this  chapter  we  describe  the  basic  idea  of  the  block  Lanczos 
method,  develop  the  algorithm,  present  several  simulation  results,  and  compare 
them  with  those  of  the  single  vector  algorithm. 


29 


A.    BLOCK  LANCZOS  METHOD 

Given  a  large  real  correlation  matrix  R  of  size  n*n,  we  can  generate  a  banded 
tridiagonal  matrix  Ts  of  size  qs*qs,  where  qs<n,  q  is  the  block  size,  and  s  is  the 
number  of  the  Lanczos  blocks.  Starting  from  an  initial  n*q  orthonormal  matrix  Qt, 
the  purpose  is  to  compute  a  sequence  of  mutually  orthonormal  n*q  matrices  Q2,  Q3, 
Q4,  •  •  • ,  Qs  such  that  the  space  of  vectors  spanned  by  the  columns  of  these  matrices 
contains  the  columns  of  the  matrices  Q^  RQl5  R2Ql5  ...,  Rs_1Ql5  where  0<q<f  and 
l<s<^.  Note  that  usually  we  have  ps<<n.  The  block  Lanczos  algorithm  can  be 
summarized  as  follows  [Ref.  1]:   For  i=l,2, •  •  -,s,  compute 


Aj  =  QT(RQ.  -  Q..^)  (3-1) 


Pi  =  RQi-QiAi-Qi.1BT  (3-2) 


Q.+1Bi+1  =  Pj     (QR  factorization  of  Pj).  (3-3) 

In  this  procedure,  Qi+1  is  orthonormal  to  all  previous  Qj.  The  purpose  of  the  block 
Lanczos  algorithm  is  to  find  Aj,  Bi+1  and  Qi+1,  where  Aj  and  Bi+1  are  the  element 
matrices  of  the  desired  tridiagonal  matrix,  and  Qi+1  is  the  (i+l)st  orthonormal 
Lanczos  block.  The  blocks  Qi:  for  i=l,2,  •  •  •  ,s,  form  an  orthonormal  basis  of  the 
Krylov  subspace,  defined  as 

JTS(Q15R)  =  spantQ^RQ^Qi,  ■  •  ■ ,  Rs_1Qi}  (3-4) 

corresponding  to  the  first  block  Qv    The  matrix  Q,    defined  as  a  catenation  of  the 
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Lanczos  blocks 


Q  =  [Qi  Q2  Q3  •  •  •  QJ 


(3-5) 


is  orthogonal  [Ref.  4],    i.e.,    Q  Q  =  I  or  Q    =  Q"1.    We  then  generate  the  banded 
tridiagonal  matrix  Ts  using  Aj  and  Bj: 


T   = 


A,    B* 
B2    A2 


B: 


B3    A3 


0 


Bs    As 


(3-6) 


The  elements  Aj  and  Bj  are  q*q  coefficient  matrices,  so  that  Ts  is  a  qs*qs  banded 
matrix.  The  off-diagonal  blocks,  Bj,  are  upper  triangular  matrices  and  the  main 
diagonal  blocks,  Aj,  are  symmetric  matrices  so  that  Ts  itself  is  a  symmetric  matrix. 
Also,  the  banded  tridiagonal  matrix  Ts  can  be  determined  by 


T   ^Q'RQ. 


(3-7) 


There  are  basically  two  different  types  of  block  Lanczos  procedures,  namely,  the 
iterative  procedure  and  the  noniterative  procedure.  The  noniterative  procedure 
follows  along  the  lines  of  the  single  vector  Lanczos  procedure.  In  this  method  a 
sequence  of  blocks  {Q1,Q2,-  •  •}  are  generated  where  the  length  of  the  sequence  can 
be  determined  by  the  size  of  Ts  and  the  amount  of  storage  available.  The  generated 
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Lanczos  blocks  may  or  may  not  be  reorthogonalized.  The  iterative  procedure  uses 
the  block  recursion  to  generate  the  block  tridiagonal  matrix.  First,  the  relevant 
eigenvalues  and  eigenvectors  of  these  block  tridiagonal  matrices  are  computed. 
Then  the  corresponding  Ritz  vectors  are  computed  and  used  as  updated 
approximations  to  the  desired  eigenvectors.  If  convergence  has  not  been  achieved  in 
k  iterations  (see  Step  2  below),  more  iterations  are  performed  these  updated 
eigenvectors  until  the  procedure  has  converged. 

The  following  steps  show  the  basic  iterative  block  Lanczos  procedure: 

Step  1.  For  k=l  start  with  an  initial  arbitrary  n*q  block  Qk  where  the 
columns  of  Qk  are  orthonormal. 

Step  2.  Compute  Pk  =  RQk  -  QkAk  using  Qk  where  Ak  =  (Qk)TRQk  and 
use  the  norms  of  the  columns  of  Pk  to  check  for  convergence.  If 
convergence  has  occured,  then  stop;  otherwise,  go  to  Step  3. 

Step  3.  Generate  a  sequence  of  blocks  Qk  using  the  recursion  in  Eqn(3— 2) 
and  Eqn(3— 3)  for  j=2,  3,  •••,  s.  Use  the  coefficient  matrices  Ak 
and  Bj+1  to  define  the  real  symmetric  block  tridiagonal  matrix  Tk. 

Step  4.  Compute  the  q  algebraically  smallest  eigenvalues  of  Tg  and  the 
correspoding  eigenvectors  Yk  where  Yk  =  {yk,  yk,  •  •  • ,  yk}. 

k  +  l 

Step  5.        Obtain  the  new  Lanczos  block  Qj     given  by 


k  +  l  k     k 

Q,     =  QYk  (3-8) 


k  k        k  k 

where  Q    =  {Qj,  Q2, •  •  • ,  Qs}.  Increment  k  to  k+l  and  go  to  Step  2. 
From  the  above  procedure  we  can  generate  the  block  tridiagonal  matrix  Tg  and  also 
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compute  the  eigenvectors  Yk  which  approximate  the  eigenvectors  of  R.    As  noted, 
Step  2  provides  a  check  for  the  convergence  of  the  procedure. 

If  it  were  the  case  that  qs  =  n,  then  Tg  would  be  similar  to  R  and  the 
eigenvalues  of  Tg  would  also  be  similar  to  the  eigenvalues  of  R.  Particularly,  some 
of  the  extreme  eigenvalues  of  Ts  would  be  approximately  the  same  as  the 
corresponding  eigenvalues  of  R.  Generally,  because  of  the  numerical  properties  of 
the  block  Lanczos  algorithm,  it  is  not  practical  to  carry  the  method  through  to 
completion  [Ref.  1].  The  importance  of  the  algorithm  lies  in  the  fact  that  some  of 
the  smallest  (and  largest)  eigenvalues  of  Ts  will  closely  approximate  the 
corresponding  eigenvalues  of  R  for  values  of  s  such  that  qs  <<  n.  This  is  stated  by 
the  following  theorem. 

Theorem  3.1  [Ref.  1].  Let  Xl  <  A2  <  •  •  •  <  An  be  the  eigenvalues  of  R  and 
vv  v2,  ••-,  vn  be  the  corresponding  orthonormal  eigenvectors.  Assume  that 
Aq<Aq+1.  Apply  the  block  Lanczos  recursion  in  Eqns(3— 1)— (3— 3)  to  R  generating  s 
blocks  and  let  //j  <  //2  <  •  •  •  <  Mn  be  the  eigenvalues  of  Ts.  Suppose  that 


W  = 


W, 

w2 


=  VTQ,  (3-9) 


is  the  n*q  matrix  of  projections  of  the  starting  block  of  vectors  Qj  on  the 
eigenvectors  of  R,  where  V  =  [vt,  v2,  v3,  •  •  •  ,  vj  and  Wj  is  a  q*q  matrix  composed 
of  the  first  q  rows  of  W.  Suppose  further  that  Wj  is  nonsingular  so  that  amin,  the 
smallest  singular  value  of  Wj,  is  greater  than  zero.  Then  for  k  =  1,  2,  3,  ...,  q,  the 
eigenvalues  of  Ts  satisfy 

Ak  <  //k  <  Ak  +  el  (3-10) 
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where  the  spread  ek  is  given  by 


(An  -Ak)tan20 


4=     — r^^  (3-H) 


l"7k 


and  $  =  cos-10-min,  7k  =  (Ak  -  Aq+1)/(Ak  -  An),  and  ^s.t  is  the  (s-l)th    Chebyshev 
polynomial  of  the  first  kind. 

This  theorem  illustrates  the  importance  of  the  local  gaps  |  Ak  —  Aq+1| ,  but 
does  not  show  the  potential  positive  effect  of  the  outer  loop  iteration  of  an  iterative 
block  Lanczos  procedure  on  reducing  the  overall  effective  spread  and  thereby 
improving  the  convergence  rate. 

Note  that,  as  we  have  defined  it,  the  block  Lanczos  method  is  not  a  method  for 
finding  the  eigenvalues  and  the  eigenvectors  of  a  symmetric  matrix  R.  Rather,  it  is 
a  procedure  for  finding  a  block  tridiagonal  matrix  Ts  which  is  similar  to  R.  To 
produce  a  complete  algorithm  for  finding  the  eigenvalues  and  the  eigenvectors,  we 
need  to  combine  the  Lanczos  algorithm  with  a  technique  for  finding  the  eigenvalues 
fi   and  the  eigenvectors  yk  of  Ts  such  as  the  QR  algorithm. 

Now,  we  will  consider  certain  properties  of  the  block  Lanczos  algorithm  and 
problems  associated  with  its  implementation  and  application. 

The  computed  Lanczos  blocks  Qj  are  desired  to  be  mutually  orthogonal.  In 
practice,  however,  because  of  the  arithmetic  errors  when  Pj  is  computed,  they 
rapidly  lose  orthogonality.  Thus,  after  a  few  iterations  of  the  block  Lanczos 
algorithm,  the  current  Qj  is  no  longer  orthogonal  to  the  previous  Lanczos  blocks 
Qi>  Q2'  ' '  '■>  Qi-i-   The  subsequent  losses  in  orthogonality  between  the  blocks  caused 
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by  the  roundoff  errors  increase  as  we  increase  the  number  of  Lanczos  blocks.  At 
some  stage  when  these  losses  have  accumulated  sufficiently,  the  assumption  that 
the  block  tridiagonal  Lanczos  matrix  is  the  projection  of  the  given  matrix  R  on  the 
subspaces  Q  will  be  false  [Ref.  1].  As  a  result  the  q  smallest  eigenvalues  of  Ts  may 
not  approximate  the  q  smallest  eigenvalues  of  R.  Thus,  it  requires  costly 
reorthogonalization  of  each  Qi+1  with  respect  to  all  the  previous  Lanczos  blocks  to 
maintain  the  stability  of  the  algorithm  [Ref.  4]. 

Loss  of  orthogonality  goes  hand-in-hand  with  the  convergence  of  some  of  the 
eigenvalues  of  Ts  to  the  eigenvalues  of  R.  In  this  case  we  have  two  options,  stop 
the  Lanczos  iterations  as  the  blocks  begin  to  lose  their  mutual  orthogonality  or 
reorthogonalize  the  blocks  if  more  iterations  are  desired.  The  difficulty  in  using  the 
Lanczos  method  in  this  way  lies  in  reliably  and  efficiently  determining  at  what  point 
the  orthogonality  is  being  lost.  In  order  to  reorthogonalize  all  Lanczos  blocks,  we 
first  reorthogonalize  the  residual  matrix  Pi  with  respect  to  all  the  previous  Lanczos 
blocks  and  compute  the  next  block  Qi+1  and  the  coefficient  matrix  Bi+1  such  that 
P~Qi+1Bi+1.  This  modification  preserves  the  stability  of  the  algorithm  but  at  a 
considerable  cost  because  the  reorthogonalization  process  requires  a  large  number  of 
arithmetic  operations.  The  need  for  reorthogonalization  seems  to  increase  with  the 
size  n  of  R  and  the  number  of  blocks  that  are  required  to  be  computed  in  the 
algorithm  [Ref.  1]. 

We  now  describe  how  to  choose  the  block  size  q.  It  is  usually  best  to  choose  q 
equal  to  the  number  of  eigenvalues  and  eigenvectors  r  that  we  are  attempting  to 
compute  which  could  be  the  multiplicity  of  the  smallest  eigenvalue.  Theorem  3.1 
suggests  that  a  good  choice  for  q  is  one  for  which  the  gap  between  Aq  and  Aq+1  is 
fairly  large. 
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B.    ALGORITHM 

It  is  possible  to  find  a  few  of  the  extreme  eigenvalues  and  corresponding 
eigenvectors  of  a  real  symmetric  matrix  using  the  block  Lanczos  algorithm  rather 
than  computing  the  entire  matrix  decomposition.  Each  of  these  smallest  eigenvalues 
and  the  corresponding  eigenvector  of  the  autocorrelation  matrix  for  received  signals 
from  a  sensor  array  has  the  spectral  information  to  estimate  the 
direction-of-arrival. 

The  computer  simulation  experimental  set  up  used  is  the  same  as  the  one  shown 
in  Figure  1.  We  receive  the  signals  at  a  linear  array  containing  M  equally -spaced 
sensors  and  generate  the  n*n  autocorrelation  matrix  of  these  received  signals. 

1.       Reduction 

The  reduction  of  the  data  proceeds  as  follows.  Using  the  block  Lanczos 
algorithm  we  reduce  the  autocorrelation  matrix  R  into  a  block  tridiagonal  matrix 
that  has  the  same  extreme  eigenvalues  as  R.  In  this  section  we  will  present  a 
practical  algorithm  to  implement  the  block  Lanczos  method. 

Given  an  n*n  autocorrelation  matrix  R,  we  generate  an  initial  n*q  matrix 
Qj  which  is  chosen  arbitrarily  and  orthonormalized.  The  number  of  vectors  in  each 
block,  q,  is  considered  to  be  between  3  and  5  in  this  study.  To  begin,  we  compute 
RQj  and  a  residual  matrix  Pj  given  by 

P,  =  RQ,  -  QjA,  (3-15) 

where  A1  is  a  q*q  coefficient  matrix  chosen  so  that  the  Euclidean  norm  of  P2  is 
minimized  with  respect  to  all  possible  choices  of  A{  [Ref.  4].  It  can  be  shown  [Ref.  1] 
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that  IIPJI  is  minimized  when 


Ax  =  qIrQj.  (3-16) 


The  qxq  matrix  Ax  forms  the  first   block  on  the  main -diagonal   of  the  block 
tridiagonal  matrix,  Ts  (see  Eqn(3— 6)).  With  this  choice  for  Av  we  have 


Pi  =  (I "  QiQT)RQi  (3-17) 


That  is,  Pt  is  the  projection  of  RQj  onto  the  space  orthogonal  to  that  spanned  by 
the  columns  of  Qv  The  second  Lanczos  block  of  vectors  Q2  and  a  qxq  upper 
triangular  coefficient  matrix  B2  are  then  obtained  by  using  the  QR  factorization 
with  modified  Gram-Schmidt  procedure  on  Px: 

Q2B2  =  Pv  (3-18) 

The  current  Lanczos  block  Q2  is  orthonormal  to  the  previous  block  Qv  The  upper 
triangular  matrix  B2  and  its  transposed  version  B2  form  the  first  elements  in  the 
sub-diagonal  and  the  super-diagonal,  respectively,  of  the  block  tridiagonal  matrix 
T 

The  remaining  matrices  in  the  sequence  of  the  Lanczos  blocks  Ql5  Q2,  •  •  •, 
Qs,  where  s<<n,   are  computed  as  follows:   For  i=2,3,-  •  -,s,  compute 


As  =  QT(RQ.  -  Q-.^T)  (3-19) 
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P^RQi-QiAi-Q^  (3-20) 


Qi+iBi+1  =  Pi  (3-21) 

where  Qi+1  and  Bi+1  are  obtained  as  the  QR  factorization  of  the  residual  matrix  Pv 
A  modified  Gram-Schmidt  procedure  can  be  used  to  reorthgonalize  the  columns  of 
Pj.  This  means  that  Qi+1  is  orthonormal  to  all  previous  matrices  Q1?  Q2,  •  •  •,  Qj. 
As  we  increase  the  number  of  the  Lanczos  blocks,  their  mutual  orthogonality  is 
preserved  because  of  the  built-in  QR  algorithm  to  factor  the  residual  matrix  Pj. 
Consequently,  the  space  spanned  by  Ql5  Q2,  •  •  •,  Qs  contains  the  columns  of  the 
matrices  Q1?  RQP  R2Qj,  •  •  •,  RS-1Q!  which  form  the  orthonormal  basis  of  the 
Krylov  subspace.  Thus,  we  do  not  need  to  reorthogonalize  the  Lanczos  blocks  when 
we  use  the  algorithm  in  Eqns(3— 17)— (3— 19). 

2.       Eigendecomposition 

We  need  to  compute  the  eigenvalues  of  the  block  tridiagonal  matrix  Ts 
which  approximate  the  eigenvalues  of  an  autocorrelation  matrix.  Then  the 
eigenvectors  of  R  corresponding  to  these  eigenvalues  are  determined  by  knowing  the 
matrix  of  Lanczos  blocks. 

There  are  several  techniques  for  computing  the  eigenvalues  and 
eigenvectors  of  a  given  matrix.  The  fundamental  algebraic  eigenproblem  is  to 
determine  the  eigenvalues  n^  given  the  set  of  qs  homogeneous  linear  equations  in  qs 
unknowns  [Ref.  17] 

VIs  Wi        fori=l,2,-.-,qs  (3-22) 
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where  the  eigenvectors  yi  of  Ts  satisfy 

0,  ifi*  j, 


T 

yjyj  = 


(3-23) 
.  1,  ifi  =  j, 


that  is,  yj  are  mutually  orthonormal  vectors.  From  Eqn(3— 22)  the  characteristic 
equation  associated  with  the  matrix  Ts  is  given  by 

det(Ts-//iI)  =  0.  (3-24) 

Expanding  the  determinant,  we  have  the  polynomial  equation 

a0  +  ajt  +  •  •  •  +  aqs^qs-1  +  //qs  =  0  (3-25) 

where  a:  are  the  coefficients  of  $  and  the  roots  of  this  polynomial  give  us  the 
eigenvalues  //j  of  the  Ts.  Corresponding  to  any  eigenvalue  //i5  Eqn(3— 22)  has  at 
least  one  non -trivial  solution  yj.  Since  the  eigenvalues  of  Tg  approximate  a  few  of 
the  extreme  eigenvalues  of  autocorrelation  matrix  R,  we  choose  the  q  smallest 
eigenvalues  and  the  corresponding  eigenvectors  of  Ts  and  use  the  Ritz  vector  to 
compute  the  eigenvectors  of  R  given  by  [Ref.  1] 

X  =  QY  (3-26) 

where  X  is  a  group  of  approximate  eigenvectors  of  R,  X=[xj  Xj  •  •  •  xq],  Q  is  a 
matrix  of  Lanczos  blocks,  and  Y  is  a  group  of  q  smallest  eigenvectors  of  Tg,  Y=[yj 
y2  '  "yj-   Also,  we  can  use  the  Rayleigh  quotient  iteration  to  find  the  eigenvectors 
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of  R  related  to  the  eigenvalues  of  Tg.    However,    this  method  is  slower  than  using 
the  Ritz  vector  because  it  needs  more  time  to  converge  for  a  given  eigenvalue. 

Now,  using  the  eigenvectors  of  the  autocorrelation  matrix  R  we  estimate 
the  spatial  spectrum  of  the  received  signals.  Each  of  the  eigenvectors  contains  the 
true  spectral  information  as  well  as  some  spurious  peaks.  The  direction -of-arrival 
of  point  sources  can  be  estimated  by  computing  the  spatial  power  spectral  density  of 
the  eigenvectors  of  R.  The  power  spectral  density  estimate  for  the  fh  eigenvector 
corresponding  to  one  of  the  smallest  eigenvalues  of  R  is  computed  as  follows: 


sj*(/)  = 


n-1 


V     x    z"1 
1=0 


**P*f 


(3-27) 


th 


where  x-j  are  the  elements  of  the  j tn  eigenvector,  X:  and  0  <  /<  0.5. 


C.    SIMULATION  RESULTS 

Using  the  block  Lanczos  algorithm  we  can  selectively  compute  a  few  of  the 
smallest  eigenvalues  and  eigenvectors  of  an  autocorrelation  matrix.  These 
eigenpairs  are  in  the  noise  subspace  and  contain  the  spectral  information  of  the 
source  bearings  from  an  array  of  sensors.  Thus,  we  could  estimate  the  spatial 
power  spectral  density  for  each  eigenvector  using  Eqn(3— 27).  We  have  used  the 
spectral  product  of  several  individual  PSDs  to  improve  the  direction-of-arrival 
estimation  performance.  The  advantage  in  using  the  multiplicative  PSD  function 
was  demonstrated  in  Chapter  2  for  the  single  vector  Lanczos  algorithm.  The  same 
advantage  holds  for  the  block  Lanczos  algorithm. 
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The  computer  simulation  experiment  consists  of  an  equally -spaced  linear  array 
of  100  sensors  receiving  signals  of  known  temporal  frequency  from  various  bearings. 
The  size  of  the  autocorrelation  matrix  computed  is  25x25.  The  number  of  smallest 
eigenpairs  to  be  estimated  is  q,  which  is  the  size  of  the  Lanczos  block.  We  choose  q 
to  be  3,  5,  or  7  in  these  examples  and  used  5  dB,  0  dB,  -5  dB  and  -10  dB 
signal-to-noise  ratios  (SNR). 

Example  1  is  the  detection  of  a  target  at  9°  for  different  SNRs.  Figure  10  shows 
the  spectral  overlay  of  5  eigenvectors  corresponding  to  the  smallest  eigenvalues  at 
an  SNR  of  5  dB.  The  spectrum  of  each  eigenvector  has  a  common  peak  at  the  true 
bearing  9°.  Figure  11(a)  and  Figure  11(b)  show  improved  DO  A  estimation  where 
the  former  figure  illustrates  the  average  of  5  eigenvectors  and  the  latter  figure  shows 
the  spectral  multiplication  of  those  eigenvectors.  As  can  be  seen,  the  spectral 
multiplication  technique  has  greatly  improved  the  nulls  and  the  peak  at  9°.  Thus, 
in  the  remainder  of  the  results,  we  have  used  the  spectral  multiplication  method  to 
compute  the  PSD  function.  Figure  12(a)  shows  the  performance  at  0  dB  using  3 
eigenvectors  and  Figure  12(b)  shows  the  result  at  an  SNR  of  0  dB  using  5 
eigenvectors.  Improved  performance  is  obtained  by  using  more  eigenvectors.  In 
Figure  12(a)  we  have  several  large  spurious  peaks  around  the  true  peak  at  9°,  but 
in  Figure  12(b)  we  have  just  2  small  spurious  peaks  and  the  true  bearing  is  clearly 
evident.  Figure  13(a)  and  Figure  13(b)  show  the  results  at  an  SNR  of  -5  dB  using  3 
and  5  eigenvectors,  respectively.  The  spurious  peaks  have  increased  magnitudes 
but  the  difference  between  the  true  peak  and  the  spurious  peaks  is  still  large  enough 
to  determine  the  true  bearing.  Figure  14(a)  and  Figure  14(b)  show  the  result  at  an 
SNR  of  —10  dB.  Using  3  eigenvectors  (Figure  14(a)),  it  is  hard  to  determine  the 
true  bearing.     Using  5  eigenvectors  (Figure  14(b)),  the  true  peak  appears  to  be 
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larger  than  the  spurious  peaks  and  we  can  recognize  the  true  bearing. 

In  Example  2  we  have  3  targets  at  34°,  36°  and  54°.  Two  targets  are  very 
closely  spaced  in  bearing.  At  0  dB,  the  results  of  two  cases,  where  one  has  used  5 
eigenvectors  (Figure  15(a))  and  the  other  7  eigenvectors  (Figure  15(b)),  indicate 
good  performance.  There  are  no  spurious  peaks  and  the  resolution  is  clearly 
acheived.  Figure  16(a)  shows  the  result  at  —5  dB  using  5  eigenvectors.  In  this 
result,  the  resolution  is  acheived  but  one  spurious  peak  is  as  large  as  the  true  peaks. 
In  Figure  16(b),  the  result  shows  improved  performance  using  7  eigenvectors. 
There  are  no  spurious  peaks  and  good  spectral  resolution  is  acheived.  At  —10  dB, 
the  algorithm  cannot  separate  the  two  targets  located  at  34°  and  36°  when  we  use  5 
eigenvectors  (Figure  17(a)).  When  we  use  7  eigenvectors  however  (Figure  17(b))  it 
can  separate  the  two  closely  located  sources.  Nevertheless,  a  number  of  spurious 
peaks  are  larger  in  magnitude  than  the  peak  at  36°  making  it  impossible  to 
determine  the  true  bearings  accurately. 

The  results  in  this  chapter  show  that  the  eigenvectors  computed  using  the  block 
Lanczos  algorithm  can  be  used  to  determine  the  spectrum  even  in  very  low  SNRs. 
Since  the  block  Lanczos  method  can  compute  a  few  of  the  extreme  eigenpairs  of  a 
large  symmetric  matrix,  it  is  efficient  in  computations  and  storage. 
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Figure  10.       One  target  at  9°  (5  dB)  :  overlay  of  5  eigenvectors 
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Figure  11.       One  target  at  9°  (5  dB)  :     (a)  average  of  5  eigenvectors 
(b)  product  of  the  PSDs  5  eigenvectors 
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Figure  13.       One  target  at  9°  (-5  dB)  :  product  of  the  PSDs 
(a)  3  eigenvectors       (b)  5  eigenvectors 
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Figure  14.       One  target  at  9°  (-10  (IB)  :  product  of  the  PSDs 
(a)  3  eigenvectors       (b)  5  eigenvectors 


47 


10° 
10-' 
10   z 
2?     10-:' 

•r-« 
W 

I    to- 

re 

S   io-5 


hi 


10 


o 
Ou      10" 


1(1    " 


10- 


I0-io 


0 


.L.  i IJ . ._. 

20  10  CO  00  100 

Hearings 


10° 


10' 

- 

1              ■  ■  ■■  - 

1                            ,......- 

JO2 

- 

- 

.~    10  3 

c 

Q    1U4 

- 

1 

- 

13 

S      105 

ex 
en     106 

- 

- 

o     107 

- 

10  R 

- 

- 

109 

- 

- 

1010 

, 

,_ 

1                             I 

20 


40  60 

Bearings 


80  100 


Figure  15.       Three  targets  at  34°,  36°  and  54°  (OdB)  : 

product  of  the  PSDs  (a)  5  eigenvectors     (b)  7  eigenvectors 
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Three  targets  at  34°,  36°  and  54°  (-5  dB)  : 

product  of  the  PSDs  (a)  5  eigenvectors     (b)  7  eigenvectors 
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Figure  17.       Three  targets  at  34°,  36°  and  54°  (-10  dli)  : 
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D.    COMPARISON 

In  this  section  we  compare  the  performance  of  the  block  Lanczos  algorithm  and 
the  single  vector  Lanczos  algorithm  for  DOA  estimation.  In  both  cases,  we  chose 
the  seven  smallest  eigenvalues  and  corresponding  eigenvectors  and  use  the  spectral 
multiplication  method  since  this  gave  the  best  preformance  for  both  algorithms. 

Example  1  is  the  comparison  for  detection  of  two  targets  at  18°  and  45°  for 
different  SNRs.  Figure  18(a)  shows  the  result  using  the  single  vector  Lanczos 
algorithm  and  Figure  18(b)  indicates  the  result  using  the  block  Lanczos  algorithm  at 
an  SNR  of  0  dB.  The  result  of  the  single  vector  case  has  two  small  spurious  peaks; 
the  block  case  has  one  very  small  spurious  peak.  In  both  cases  the  true  bearings  are 
clearly  distinguished  from  the  spurious  peaks.  Figure  19  shows  a  comparison  of  the 
results  at  —5  dB.  The  performance  of  the  block  Lanczos  algorithm  (Figure  19(b))  is 
much  better  than  the  single  vector  case  (Figure  19(a))  even  though  several  large 
spurious  peaks  are  present.  Figure  20(a)  is  the  result  at  —10  dB  using  the  single 
vector  Lanczos  algorithm  and  Figure  20(b)  is  the  corresponding  result  of  the  block 
case.  In  the  single  vector  case,  the  spurious  peaks  are  almost  the  same  as  the  true 
peaks  and  we  cannot  distinguish  the  true  bearing  location.  In  the  block  case, 
however,  the  true  peaks  are  slightly  larger  than  the  largest  spurious  peak. 

Example  2  is  a  comparison  for  detection  of  four  targets  at  18°,  27°,  29°,  and 
45°.  Figure  21  shows  the  performance  of  the  two  methods  at  0  dB.  The  resolution 
is  clearly  achieved  in  both  cases.  Figure  22  shows  the  results  at  —5  dB. 
Figure  22(a)  shows  that  the  single  vector  Lanczos  algorithm  cannot  separate  the  two 
closly  spaced  targets.  The  block  case  (Figure  22(b))  shows  slightly  better  resolution 
between  targets  located  at  27°  and  29°.  At  -10  dB  (Figure  23(a)  and  Figure  23(b)) 
neither  method  can  separate  the  two  closely  spaced  in  bearings  at  27°  and  29°. 
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Also,  a  spurious  peak  is  produced  by  both  methods. 

The  block  Lanczos  algorithm  is  known  to  estimate  the  multiple  eigenvalues 
better  than  the  single  vector  case  [Ref.  1].  This  situation  is  applicable  to  the  DO  A 
estimation  using  noise  subspace  computation  where  the  noise  is  white.  As  observed 
in  the  results  of  Figures  18  —  23,  the  spectral  estimation  performance  of  the  block 
Lanczos  algorithm  is  consistently  better  than  the  single  vector  algorithm, 
particulary  at  low  SNRs.  Also,  the  block  method  provided  better  spectral 
resolution  than  the  single  vector  method  in  our  tests.  While  more  analysis  is  needed 
to  validate  the  simulation  results,  the  overall  performance  of  the  block  method  is 
quite  encouraging. 
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Figure  18.     Two  targets  at  18°  and  45°  (0  (IB)  : 

(a)  single  vector  Lanczos  method     (b)  block  Lanczos  method 
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Figure  19.     Two  targets  at  18°  and  45°  (-5  dI3)  : 

(a)  single  vector  Lanczos  method     (b)  block  Lanczos  method 
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Figure  20.     Two  targets  at  18°  and  45°  (-10  (IB)  : 

(a)  single  vector  Lanczos  method     (b)  block  Lanczos  method 
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Four  targets  at  18°,  27°,  29°  and  45°  (0  dH)  : 

(a)  single  vector  Lanczos  method     (b)  block  Lanczos  method 
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Figure  22.     Four  targets  at  18°,  27°,  29°  and  45°  (-5  (IB)  : 

(a)  single  vector  Lanczos  method     (b)  block  Lanczos  method 
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IV.    SUMMARY  AND  CONCLUSIONS 

In  this  thesis  we  examined  the  use  of  the  single  vector  Lanczos  method  and  the 
block  Lanczos  method  and  its  application  to  spectral  analysis  and 
direction -of -arrival  problems. 

We  computed  a  few  of  the  extreme  eigenvalues  and  their  associated 
eigenvectors  of  a  large  symmetric  matrix  using  the  block  Lanczos  method.  The 
eigenvalues  and  eigenvectors  of  the  Lanczos  matrix  Tg  approximate  the 
corresponding  eigenvalues  and  eigenvectors  of  the  given  matrix  R.  The  block 
Lanczos  algorithm  can  directly  determine  the  multiplicities  of  the  effective 
eigenvalues  and  the  eigenvectors  of  R.  We  found  that  the  spectral  estimate  of  the 
block  Lanczos  method  is  more  accurate  than  the  single  vector  Lanczos  method, 
particulary  at  low  SNRs.  Since  we  compute  only  a  few  of  the  extreme  eigenpairs  of 
a  large  autocorrelation  matrix,  the  result  of  this  algorithm  is  savings  in 
computations  and  storage.  This  algorithm  may  be  applied  to  any  system  where  one 
needs  to  obtain  a  rapid  decomposition  of  a  large  correlation  matrix. 

Although  the  results  of  this  thesis  are  most  encouraging,  some  additional 
work  still  remains  to  be  done.  We  need  to  compare  the  results  in  computational 
speed  and  accuracy  with  other  eigendecomposition  techniques  for  valdating  this 
algorithm.  Also,  we  need  to  analyze  the  effect  of  roundoff  errors  on  the  eigenpair 
estimation  of  the  Lanczos  algorithm. 
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